package ru.usu.gv.utils.math.mds;

import ru.usu.gv.utils.math.Matrix;
import ru.usu.gv.utils.math.Vector;

import java.util.ArrayList;
import java.util.List;

/**
 * @author spupyrev
 * 22.11.2009
 * 
 * Computes the mds by distance scaling
 */
public class DistanceScalingTest
{
	private static final int MAX_ITERATIONS = 10000;
	private static final double EPS = 1e-10;

	public int mds(Matrix d, Matrix w, int dimensions)
	{
		return mds(d, w, dimensions, true);
	}

	public int mds(Matrix d, Matrix w, int dimensions, boolean useClassicalInitialization)
	{
		int n = d.getRowCount();
		Matrix res;
		if (useClassicalInitialization)
		{
			res = new ClassicalScaling().mds(d, dimensions);
		}
		else
		{
			res = new Matrix(n, dimensions);
			for (int i = 0; i < dimensions; i++)
			{
				Vector v = new Vector(n);
				v.initRandomUnitLengthVector(n, i);
				for (int j = 0; j < n; j++)
					res.set(j, i, v.get(j));
			}
		}

		double[] wSum = new double[n];
		int[][] wNN = new int[n][];
		for (int i = 0; i < n; i++)
		{
			List<Integer> notNullIndexes = new ArrayList<Integer>();
			for (int j = 0; j < n; j++)
			{
				if (i != j)
				{
					if (w.get(i, j) != 0.0)
					{
						wSum[i] += w.get(i, j);
						notNullIndexes.add(j);
					}
				}
			}

			if (wSum[i] == 0)
				System.err.println(i);
			assert (wSum[i] != 0);

			wNN[i] = new int[notNullIndexes.size()];
			for (int j = 0; j < notNullIndexes.size(); j++)
				wNN[i][j] = notNullIndexes.get(j);
		}

		int iter;
		for (iter = 0; iter < MAX_ITERATIONS; iter++)
		{
			boolean changed = false;
			for (int i = 0; i < n; i++)
			{
				Vector resNew = new Vector(dimensions);
				for (int t = 0; t < wNN[i].length; t++)
				{
					int j = wNN[i][t];
					double inv = 0;
					for (int k = 0; k < dimensions; k++)
					{
						double value = res.get(i, k) - res.get(j, k);
						inv += value * value;
					}

					inv = Math.sqrt(inv);
					if (inv > EPS)
					{
						inv = 1.0 / inv;
						for (int k = 0; k < dimensions; k++)
						{
							double value = w.get(i, j) * (res.get(j, k) + d.get(i, j) * (res.get(i, k) - res.get(j, k)) * inv);
							resNew.add(k, value);
						}
					}
				}

				for (int k = 0; k < dimensions; k++)
					resNew.set(k, resNew.get(k) / wSum[i]);

				for (int k = 0; k < dimensions; k++)
				{
					if (Math.abs(resNew.get(k) - res.get(i, k)) > 1e-1)
					{
						changed = true;
					}
					
					if (resNew.get(k) != res.get(i, k))
					{
						res.set(i, k, resNew.get(k));
						//changed = true;
					}

				}
			}

			if (!changed)
				break;
		}

		//System.out.println("Stress: " + stress(d, w, res));
		return iter;
	}

	private double stress(Matrix d, Matrix w, Matrix res)
	{
		int n = d.getRowCount();
		double sum = 0;
		for (int i = 0; i < n; i++)
			for (int j = 0; j < n; j++)
			{
				double dist = 0;
				for (int k = 0; k < res.getColumnCount(); k++)
				{
					double value = res.get(i, k) - res.get(j, k);
					dist += value * value;
				}

				dist = Math.sqrt(dist);

				double diff = d.get(i, j) - dist;
				sum += w.get(i, j) * diff * diff;

				//sum += w.get(i, j) * (-2.0 * dist + dist * dist);
			}

		return sum;
	}
}
